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The study of free piston Stirling engine (FPSE) requires both accurate thermodynamic and dynamic mod¬ 
elling to predict its performances. The steady state behaviour of the engine partly relies on non linear dis¬ 
sipative phenomena such as pressure drop loss within heat exchangers which is dependant on the 
temperature within the associated components. An analytical thermodynamic model which encompasses 
the effectiveness and the flaws of the heat exchangers and the regenerator has been previously developed 
and validated. A semi-analytical dynamic model of FPSE is developed and presented in this paper. The 
thermodynamic model is used to define the thermal variables that are used in the dynamic model which 
evaluates the kinematic results. Thus, a coupled iterative strategy has been used to perform a global sim¬ 
ulation. The global modelling approach has been validated using the experimental data available from the 
NASA RE-1000 Stirling engine prototype. The resulting coupled thermodynamic-dynamic model using a 
standardized description of the engine allows efficient and realistic preliminary design of FPSE. 

© 2010 Elsevier Ltd. All rights reserved. 


1. Introduction 

One of the last and promising evolutions of the Stirling engine is 
the free piston mechanical arrangement designed by Beale at the 
end of the 60st [1]. The main advantages are known as, a simpler 
mechanical design, no lateral loads which reduces wear and a fol¬ 
lowing extended lifetime compared with classical Stirling engines. 
Therefore, applications such as radioisotope generator for deep 
space mission [2,3] are currently under development. 

Beside its advantages, the optimization of FPSE is a difficult task. 
Indeed, moving elements are driven by both the working fluid and 
gas springs pressures. The strokes and phase angle are then set by 
the coupled effect of the dynamic and thermodynamic parameters. 
As the volumes of the chambers are modified with the displace¬ 
ments of the piston and displacer, so are the pressure and the pres¬ 
sure losses through the heat exchangers. Moreover, the coupling 
intensity is related to the temperatures of the fluid within the en¬ 
gine. Hence, thermodynamic parameters of the engine have to be 
defined prior to any dynamic model of FPSE. 

Due to the complexity of this analysis, the isothermal assump¬ 
tion is usually adopted. Furthermore, temperatures within the en¬ 
gine are given as input parameters. Linearization methods are then 
used in the vicinity of an operating point to obtain an estimation of 
the performances of the engine [4-8]. As a result, the behaviour of 
FPSE can be evaluated. Those analyses are limited to starting 
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behaviour and the steady state characteristics can not be accu¬ 
rately modelled with linear approaches. 

Ulusoy in [9] has studied the behaviour of FPSE using perturba¬ 
tion methods. Averaging or multiple scales methods have been 
used in many studies for simplification purpose of the effects of 
nonlinear phenomena which occur in many problems related to 
engineering and science [10]. The centre manifold approach is 
the one which leads to a semi-analytical simplified model. It ac¬ 
counts for the local bifurcation behaviour in the neighbourhood 
of a fixed point of the nonlinear system [11]. The centre manifold 
approach reduces the number of equations of the original system 
in order to obtain a simplified system without loosing the dynam¬ 
ics of the original system as well as the effects of nonlinear terms 
[12]. These studies have outlined the drastic effect of dissipative 
and non-linear phenomena on the FPSE behaviour. An equivalent 
global evaluation of the losses has been used and a straight relation 
between pressure loss and the exchangers characteristics i.e. 
length, hydraulic radius and free flow area can not be established. 
Consequently, these approaches appear to be unsuitable for preli¬ 
minary design purpose and restricted to post analysis of existing 
design with given temperatures distribution. Therefore, there is a 
need for an accurate global approach including a thermodynamic 
isothermal model which can be used in accordance with dynamical 
analysis of the FPSE for preliminary design purpose. 

Pressure drop loss evaluation is based on flow rate through each 
heat exchanger and regenerator. As far as isothermal assumption 
can be made, Organ in [13] has suggested a method to evaluate 
mass flow rates within the engine. In addition to geometric 
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Nomenclature 



A 

area (m 2 ) 

%xdxp 

displacer - piston phase angle (rad) 

d 

intrinsic damping 

n 

efficiency 

dh 

hydraulic diameter (m) 

Ap 

pressure loss (N m -2 ) 

d w 

fiber diameter of the matrix regenerator 

K 

swept volume ratio 

e 

regenerator efficiency (-) 

( 

expansion chamber to heat source temperatures ratio 

/ 

frequency (Hz) 

3 

cooler to heater heat exchange coefficients ratio 

If 

friction factor (-) 

T 

expansion to compression chamber temperatures ratio 

h 

convection coefficient 

r 

heat source to sink temperatures ratio 

M 

total mass of the working fluid (kg) 



m 

mass of a sliding element (k) 

Subscripts 

P 

power (W) 

bp 

piston bounce space 

P 

pressure (N m -2 ) 

bd 

displacer gas spring space 

Lx 

length of exchanger (m) 

C 

compression chamber 

FItx 

number of tubes (-) 

CC 

compression chamber dead volume 

r 

working fluid ideal gas constant 

d 

displacer 

r h 

hydraulic radius of exchanger (dp, = 4 r hi ) 

dc 

compression side of the displacer 

s 

moving element stroke (m) 

de 

expansion side of the displacer 

T 

temperature (K) 

diss 

dissipative relative term 

V 

volume (m 3 ) 

E 

expansion chamber 

Vsw 

swept volume 

if 

free flow 

X 

moving part position (m) 

HC 

expansion chamber dead volume 



h 

heater 

Greek symbols 

k 

cooler 


porosity of the matrix regenerator 

P 

piston 

P 

dynamic viscosity (N m -2 s) 

R 

regenerator 

CO 

pulsation (rad s -1 ) 

th 

thermal relative term 

a 

swept volume phase angle (rad) 

w 

wetted 


characteristics, they are linked to the operating frequency, pistons 
strokes, phase angle and chambers temperatures. 

Temperatures within the engine are related to the heat 
exchangers performances and flow rate eventually. Representative 
thermodynamic models can be obtained by given heat source and 
sink temperatures instead of chambers ones [14]. Consequently, 
associated dynamic-thermodynamic studies can be used to obtain 
realistic results for FPSE models. 

In a first part, a standardized description of FPSE for the thermo¬ 
dynamic and dynamic approaches is given. Then, the main steps of 
the proposed thermodynamic approach are recalled. 

The evaluation of the mass flow rate and the related Reynolds 
number for the exchangers is detailed. It appears to be the corner 
stone of the approach. Values for heat transfer coefficient as well as 
the friction factor are given by experimental correlations [15]. 
Thus, it is possible to establish a strong relation between thermo¬ 
dynamic-dynamic model and design variables. Then, we recall the 
basis of the Flopf bifurcation analysis used to asses the dynamic 
analysis. The results of this semi analytical analysis are presented. 
Therefore, dynamic parameters can be tuned to match the thermo¬ 
dynamic ones. Finally, a joint thermodynamic-dynamic modelling 
strategy is proposed. The model is validated by comparison with 
the RE-1000 experimental data [16]. 


2. Analysis 

2.2. Standardized description of FPSE 

For any thermodynamic analysis, expansion and compression 
swept volumes (V swE ,V swC ) as well as their phase angle (a) are sup¬ 
posed to be known. These parameters are forced by mechanical 
kinematics in classical Stirling engines. Mean pressure (p me an). 
operating frequency (co), upper and lower external temperatures 
(' T h ,T l ) are set as control parameters. Depending on the heat 


exchangers characteristics and the given working fluid, chamber 
temperatures T E and T c are to be determined. Engine performances 
can be deduced eventually. 

On the contrary, in the case of a dynamical analysis of FPSE, 
chamber temperatures (T E ,T C ) are given. Mean pressure (p mea n) as 
well as working fluid are set. The problem at stake is to determine 
the operating pulsation (co), amplitudes and phase angle of dis¬ 
placer and piston s d , s p and a xdxp respectively. Swept volumes V swEl 
V swC as well as their phase angle ot are deduced and the engine per¬ 
formances can be evaluated eventually. 

The classical Schmidt analysis [17] and dynamical model of 
FPSE which can be found in the literature use different definitions. 
The proposed preliminary design method presupposes to run both 
thermodynamic and dynamic analyses. Thus, a common set of 
parameters is required aiming at a global approach. 

Relations between swept volumes and strokes amplitudes and 
phase are to be specified. For thermodynamic analysis, we set: 

Vr(t) = Vs W e/2(1 + cos(cot)) 

V c (t) = V swC /2(l + cos(cot - a)) 

As harmonic periodic displacements are obtained for FPSE, swept 
volumes can be alternatively given. Setting, x d = s d cos( cot) and 
x p = s p cos (cot + a xdxp ), instantaneous volumes of the expansion and 
compression chambers are: 


V E (t) = A de x d (t) 

V c (t) = A p x p {t) -A dc x d (t) 


Comparing Eqs. (1) and (2), one can deduce: 

V swE — 2 A de s d 

V swC 


2 v (Ap$p cos(oc xd xp) T" A dc s d ) + ( A p s p sin(oCxdxp)) 


a 


7i + arctan 


ApSp sm((x xdxp ) 


A dc s d + ApSp cos (oc xdxp )_ 
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The total swept volume denoted V sw can be expressed from V swE , 
Vswc and a. The swept volume ratio k = V swC lV swE is usually defined. 
Thus, total swept volume is: 



V sw£ \ 1 +/c 2 + 2/ccos(a) 



Alternatively, it can be expressed as: 


P = P 




mean j_ A dc /z-A de _L y 

1 _h T £ s T c s A P 


wherein: 



From the dynamic analysis, using (3), we can establish: 


k — Ad e Sd/ y ( ApSp cos ((Xxdx P ) T Adc^d) + ( ApSp sm((x X d xp )) 

V sw — 2AdeSd 


X 

N 


1 +K 2 - 2 K 


1 


1 +A 2 p s 2 p sin(a xdxp ) 2 /(ApSp cos(a xdxp ) 


AdcSdY 


( 5 ) 


2.2. Thermodynamic analysis 

A general analysis scheme and associated parameters for any 
Stirling engine are presented in Fig. 1. In this section, main results 
of the thermodynamic study developed in [14] are given. This anal¬ 
ysis relies on the usual Schmidt analysis (framed scheme of Fig. 1) 
which allows expressing the indicated power as a function of 
chamber temperatures. In addition, a thermal model is defined in 
order to take into account the heat exchanger effectiveness as well 
as the regenerator reheat loss (see Fig. 1). It is possible to express 
the thermal power of the engine and the operating point can be 
determined by the equality between the two expressions. 


2.2.2. Results from the Schmidt analysis 

The Schmidt analysis results expressed with both the thermo¬ 
dynamic and dynamic parameters are recalled hereafter. 

The instantaneous pressure is: 


P Pmean 1 + C OS{(Dt - 6) 



P = 


V / 2T7CC0S(a) + T 2 + 7C 2 


tan(0) 


2v + K + T 

/csin(a) 
/ccos(a) + t 


T = Tc/T e 


_ ^ cc VrTc Vhc 

VswE V S wE Tr V S yvE 

K = V swC/V swE 



Consequently, the main characteristic parameters of the engine can 
be given by the following analytic expressions: 




In which M is the total mass of gas within the working spaces and r 
the working fluid ideal gas constant. 

The indicated power is expressed as: 


Co Vs W e k(t - l)sin(a) 

271 2p y 7 ! 2 + tc 2 + 2?dcos(a) 


( 10 ) 


The ideal thermodynamic efficiency is the Carnot efficiency: 

rj i = \ -T c /T e (11) 


2.2.2. Thermal approach 
• Thermal power 

The cycle average power turns out to be given in an alternative 
way. According to the power balance of the machine as described 
as in Fig. 1, thermal power can expressed as: 



Fig. 1. Scheme for the thermodynamic analysis with heat exchangers and regenerator. 
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Pth = h k A wk T H (i. + $ri - 



In which = T E /T H is the temperature ratio between the heat source 
and the highest temperature of the engine and r = T L /T H . 

From a simple convection model and given wetted area of 
the cooler and the heater (A wk and A wh respectively), 3 can be ex¬ 
pressed as: 


hk^wk 

hh^wh 



In which h h and h k are the convection heat transfers. 


• Thermal efficiency 


The thermal efficiency of the engine is defined by the ratio of 
the available power by the added thermal power: 


*1th 


1 + 3y - { - {t<S 

1 - £ - (Ccond + (1 - e)M R C v )/(h h A wh )i( 1 - t) 



wherein M R is the fluid mass rate inside the regenerator, e is the 

regenerator effectiveness, and C v the specific heat for isochoric 
evolution. 


Table 1 

Parameter values of run #1011. 


Heat source T H [°C] 

600 

Cold source T c (°C) 

25 

Phase angle a (°) 

106.44 

Swept volume ratio k (-) 

1.022 

Swept volume V sw (cm 3 ) 

74.90 

Working fluid 

Helium 

Mean pressure p mean (MPa) 

7.034 

Overall efficiency r\ tm (%) 

23.7 

Operating frequency/(Hz) 

30.1 

Output power (W) 

955 

Wire mesh regenerator 1 


Regenerator length (mm) 

64.46 

Wire diameter d w (pm) 

88.9 

Porosity of matrix (-) 

0.759 

Regenerator hydraulic radius 

0.07 

Material conductivity 

16.3 

r hR (mm) 


(WirT 1 K -1 ) 


Displacer 1 


Rod diameter (mm) 

16.657 

Displacer weight (g) 

426 

Gas spring volume (cm 3 ) 

31.79 


X 9 optim 



Fig. 3. Evolution of the operating point for various frequencies. 


Though, the temperature ratio can be an optimization criteria 
because it is related to the efficiency, both power and efficiency 
have to be considered. Fig. 4 plots the power-efficiency curve for 
different frequencies. The experimental operating frequency of 
30 FIz of the RE-1000 run #1011 reference can be clearly identified 
as the optimal efficiency peak. 


• Thermodynamic conditions 


2.3. Analysis of the heat exchangers 


From the second law of thermodynamics and an engine opera¬ 
tion conditions, we set two constraint equations: 


n th c i - t 
Pth > o 



Therefore, the permissible values of f must be situated within 
hatched domain shown on Fig. 2. 

As an example, result from experiment run #1011 of the RE- 
1000 [16] (see Table 1) is presented in Fig. 2. It is clearly located 
in the hatched region close to the Carnot efficiency limit. 

As £ ratio must be defined as a single value, we choose here an 
optimal case. The inequalities Eq. (15) can be switched to equality. 
Consequently, the optimal ratio £ optim between the heat source 
temperature and the temperature of the expansion chamber of 
the engine can be formally written as: 



(<5 + 1)t 


_ sr + T _ 

(Ccond + (1 - e)M R C v )/(h h A wh )( 1 - t) 2 



Finally, the operating point can be obtained matching the values 
of the indicated and thermal power given by (10) and (12) respec¬ 
tively. Fig. 3 plots curves of temperature ratios t and £ optim for var¬ 
ious operating frequencies/= 2nlco. 

A discrepancy of 1.8% for the t temperature ratio and 0.23% for 
the C temperature ratio can be evaluated. 


The heat exchanger analysis appears to be a crucial facet of the 
engine model. Stirling engine performances rely on the heat trans¬ 
fer effectiveness. Moreover, the fundamental role of the regenera¬ 
tor itself leads to specific studies [18,19]. The dynamic aspect of the 
exchangers with respect to non linear dissipative phenomena has 
to be taken into account in a dynamic modelling of FPSE as a key 
effect as well. 

For the thermal and dynamic points of view, usual experimental 
correlations can be used. Thus, the Colburn Factor is related to the 
heat transfer effectiveness whereas the friction factor allows pres¬ 
sure drop evaluation. These correlations are based on the Reynolds 
number which has to be determined for each of the exchanger. This 
evaluation is detailed hereafter. 

2.3.1. Flow rate and reynolds evaluation 

The mass flow rate can be estimated from de Schmidt analysis 
as detailed in Organ work [13] dedicated to the study of the regen¬ 
erator problem. Following the generic schematic of a Stirling en¬ 
gine (see framed scheme in Fig. 1), the mass of gas M x occupying 
a region between expansion piston face and any significant section 
of the engine (e.g. heat exchanger opposite end face) can be ex¬ 
press as: 

M x = l fx(V swE /2(\ + cos(fflt)) + V„) (17) 

r lc 
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P, h [kW] 



Whereas the total mass of gas in the engine can be evaluated as: 

M = y0 + P cos{cot - 9)) (18) 

Consequently, the mass fraction ratio results in a simple expression: 

M*_J _ A + Bcos(cot) 

M ~ sT c 1 + Ccos(cot) + D sin(cot) 

As an example the mass held between the displacer and the regen¬ 
erator face of the heater: 


A = t(V swE /2 + V HC + V h ) 
B = tV swE /2 
C = [3 cos 6 
D = -ft sin 6 



If there is no volume left in the compression chamber (V c = 0), the 
whole mass of fluid is contained within the remaining spaces which 
is the case for the circle of the upper dashed line of Fig. 5. When the 
schematic displacer reaches its upper position, there is no mass of 
fluid in the expansion space (V E = 0) which is the case for the circle 
of the lower dashed line. 

Note that for the RE-1000 the compression chamber dead vol¬ 
ume represents about twice the cooler volume. The schematic on 
the right of Fig. 5 is the attempted representation of the engine 
volumes. 

From the evaluation of the masses, it is possible to give a mean 
value of mass flow rate for each of the exchanger and the regener¬ 
ator. The regenerator case is detailed hereafter. 

The interactive mass for the regenerator M Rint is defined by the 
maximum amount of the fluid mass between the displacer face and 
the compression output of the regenerator and the minimum value 
between the displacer and the regenerator expansion side. Theses 
masses are underlined in Fig. 6 which represents the evolution of 
M r and M h , in spaces {V E + V HC + V h } and [V E + V HC + V h + V R ] respec¬ 
tively. The horizontal dashed line gives the value of The res¬ 

ident mass of the regenerator is approximated by its mean value as 
(dotted-dashed line in Fig. 6): 


M 


R(mean ) 


P 


mean 


J_l0g(l/T) 
T e 1 -T 




Finally, the mean mass flow rate through regenerator is: 

Mr = — (M R (int) — ) (22) 

/ c 


Which enables to evaluate a Reynolds number as: 



Mr 4 rh R 
M Aff R 


(23) 


Therefore, it is possible to give the amount of fluid in the different 
chambers with respect to the position of the piston and displacer. 
Fig. 5 plots the mass of fluid within the different parts of the engine. 


Fig. 7 plots the Reynolds number for the regenerator as a func¬ 
tion of the temperature ratio t from the RE-100 parameters for dif¬ 
ferent operating frequencies. Due to the very small hydraulic 



Fig. 5. Evolution of the mass of fluid within the different chambers of the RE-1000. 
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Fig. 7. Evolution of the Reynolds numbers for heat exchangers and regenerator 
with respect to temperature ratio for three different operating frequencies. 


radius, a laminar flow is obtained for the regenerator. For the hea¬ 
ter and the cooler the usually reported turbulent behaviour of 
mean flow is predicted as well. 

Once the Reynolds numbers are evaluated for each exchanger, it 
is possible to use the experimental correlations as a simple way to 
deduce both their heat transfers and friction behaviours. As a key- 
point for the dynamic analysis we focus here on the later. 

2.3.2. Pressure drop coefficient 

Different from de work of De Monte [4,5], evaluations of indi¬ 
vidual pressure drops related to each of the heat exchangers are 
seeking here. The friction factor can be estimated [15] from Rey¬ 
nolds numbers. Thus, it is possible to give an estimation of the 
pressure drop in heater, cooler and regenerator denoted A ph , A pk , 
and A pR respectively, with respect to their geometrical parameters 
but also on the operating parameters V SWl k , a and co. Pressure 
drops can be calculated using the friction factor coefficient related 
to the Reynolds number. As a consequence the impact of the geo¬ 
metric parameters on the engine losses and its operation and per¬ 
formances can be studied. 

Fig. 8 plots the relative pressure drop with respect to the oper¬ 
ating frequency for each heat exchanger and regenerator. 

2.4. Dynamic analysis 
2.4.1. Governing equations 

Assuming that the space coordinates x p and x d are relative to the 
midpoint of the stroke for each moving part (piston, displacer). 
Besides, the case movement is neglected. The equations that 
describe the dynamic behaviour of an FPSE can be written for each 
of the moving part such that: 


m d x d + d d x d = ( A bd p bd - A de p E + A dc p c ) 
fflpXp + {dp + d em )x p = {App E p — ApP c ) 

with the following initial conditions: 

X PU=°. x Pl fc0 =°; M,=o=°-’ = 0 (25) 

The pressures p El p c are the pressure in expansion and compression 
spaces as defined in Fig. 9. p bd , p bp are pressures in the gas spring 
and bounce space of the engine. The effect of the mechanical dissi¬ 
pative phenomena is taken into account by means of the d p and d d 
damping coefficients. 

This model would be validated using the data from the RE-1000 
experiments. A dashpot load is used as a control parameter that 
determined the piston, displacer and pressure amplitudes and 
the engine frequency. As a consequence, we choose to define a sim¬ 
ple control parameter as a viscous damping associated to the pis¬ 
ton and denoted di oad = d p + d em . 

In Eq. (24) x p ( x d ) is the position of the piston (displacer) with 
respect to its rest equilibrium position. x p ( x d ) is the derivative of 
x p ( x d ) with respect to time. 

The instantaneous pressure p within the chambers is analyti¬ 
cally expressed as in (7) which is more suitable than (6). Besides, 
pressures of the chambers are linked by the pressure losses which 
occur in the heat exchangers as well as in the regenerator. Thus, we 
set: p E = p + (A p h + A p R + A p k ) and p c = p. 

Pressure in the gas spring and the buffer space may be evalu¬ 
ated by assuming the ideal gas relation for an adiabatic process: 

(26) 


Therefore, the general form of Eq. (24) can be given in the following 
way: 


~x d ~ 

*d 

X p 

Xp. 


- 0 

-S d d 

0 

_—Spd 


1 

—Ddd 

0 

-D pd 


0 

-Sd P 

0 

-Spp 


0 

-Dd P 

1 

-D PP 



~x d ' 


- 0 - 


x d 

+ 

m d 


Xp 

0 


-Xp- 


-1 

( 
_1 



with initial conditions given by Eq. (25). 

The expressions of the pressure losses A pe , A pk , and A pR are non¬ 
linear relations in which the classical Reynolds number is the main 
variable. Pressure terms are approximated evaluating Taylor’s 
development so is the pressure expression (7). The nonlinear stiff¬ 
ness and dissipative terms are grouped in fNL d and fNL p . 


2.4.2. Solution procedure 

Eq. (27) can also be seen as two nonlinear oscillators coupled by 
mechanical forces. A classical supercritical Flopf bifurcation defines 
the steady oscillations in such coupled oscillators. At the fixed 
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Fig. 9. Scheme for the dynamic analysis of FPSE. 



Fig. 10. Evolution of real and imaginary parts of the system’s poles with respect to 
the reduced damping coefficient d p = d /0Qd /d; 0admr . 


point also called the bifurcation point denoted x 0 , a couple of com¬ 
plex conjugated eigenvalues associated to the system crosses the 
imaginary axes. Hence, this point becomes unstable and because 
of the nonlinearities of the system a stable limit cycle can occurs 
in its vicinity as can be seen in Fig. 10. The effect of the load of 
the system represented by the damping coefficient d p is reflected 
in Fig. 10 in which its increasing leads to the stabilization of the 
system. Thus, the load may have to be controlled during the start¬ 
ing phase. 

• Unstable fixed point 

The fixed point [x d ;k d ;x p ;x p ] = [0;0;0;0] is obtained by solving 
the non-linear static equations. The evaluation of the eigenvalues 
of the constant Jacobian matrix at the fixed point q 0 allows the sta¬ 
bility analysis. They are defined as roots of the fourth-degree char¬ 
acteristic polynomial such as: 


r 4 + ocr 3 + fir 2 + yr + 3 = 0 (28) 

in which: 


oc = Dpp + D dd 

P — Dpp D dd — D pd D dp + S dd + S pp (29) 

y = -D dp S pd — D pd S dp + D dd S pp + D pp S dd 

& = Spp^dd ~ $pd$dp 

Depending on the operating parameters, the only situation for 
which oscillations can occur is: two complex conjugate roots with 
positive real parts and two complex conjugate roots with negative 
real parts. Thenceforth, a further analysis is required to asses the 
steadiness of these oscillations. 

The roots of (28) evaluated at the fixed point can be expressed 
as two complex conjugates expressions: 


Ajci = 0Ci — j&h; Jjc2 = ai + joji (30) 

2/Si = oc 2 —JCO 2 ; Aj S 2 = OC 2 +jt02 (31) 

in which j is the complex number such as j = T and ai is any po¬ 
sitive real number. 


• Centre manifold approach 

The previous linear analysis allows defining the condition nec¬ 
essary to the engine operation but can not represent the behaviour 
of the engine after the starting point. Therefore, another theory has 
to be used to study the post bifurcation behaviour. This section 
summarizes the application of the centre manifold method [12] 
to study the system defined by (27). The stability of the centre 
manifold is equivalent to the stability of the original system in 
the vicinity of the fixed point. Again, the bifurcation appears when 
one or several eigenvalues cross the imaginary axis in the complex 
plane with the variation of a control parameter (see Fig. 10). 

At the Hopf bifurcation point, Re(7j c i) = 0 and Eq. (32) can be 
written in the canonical form: 
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•; Vc] c v c + G 2 {v c , v s ) + G 3 (v c , v s ) 

; V K = Jj»s + H 2 (I' c , V s ) + H 3 (|/ c , V s ) 

wherein J c et J s have eigenvalues such as X jc = (x^+l-jco^ and 
%s = a 2 +/—^ 2 . G 3 , H 2 and H 3 are vector for which the two com¬ 
ponents are polynomials of degree 2 and 3 in the components of 
vectors v c and v s namely v c \, v c2 and v sl , v s2 . 

For the FPSE case, it may be assumed that a 2 is negative. v s de¬ 
fines a stable subspace which is the linear approximation to the 
stable manifold. When we choose an initial condition on this stable 
manifold sufficiently close to the fixed point the solution curve will 
go toward the fixed point. In accordance to the parameters, an 
unstable subspace defined by the eigenvalues of matrix A such as 
their real part is zero can exist. The centre manifold theorem al¬ 
lows to represent locally the centre manifold as {[v c , v s ] such that 
v s = h(v c ), h(0) = 0, Dh(0) = 0} and consequently reduce the four 
dimensional initial problem to a two dimensional one. 
Substituting into the second line of Eq. (32) one obtains: 

■; V=) c Vc + G 2 (v c , Vs) + G 3 (v c , v s ) (33) 

D»c(h(v c ))-;v =J s h(» c ) + H 2 (i'c,h(» c )) +H 3 (i> c ,h(i/ c )) (34) 

Thus, it is possible to define an approximate solution h using a 
power expansion without constant and linear terms: 

m p p 

h(*0 = (35) 

p=i+j= 2 i=0 j= 0 

Replacing ;v c , the complex coefficients a y can be obtained. As a re¬ 
sult, the dynamic behaviour of the system is determined by the fol¬ 
lowing reduced system: 


V = ] c v c + G 2 (v c ,h{Vc)) + G 3 (v c ,h(Vc)) 



• Normal form analysis 

The normal form analysis aims at a transformation of a system 
of nonlinear equations through a sequence of nonlinear-near iden¬ 
tity transformations to eliminate as many nonlinear terms as pos¬ 
sible. Those which cannot be removed are called the secular or the 
resonant terms. This simplest form of the equations is called the 


xj and x p [mm] 


i 


0 


7T 


2 tt 


CO t 



“normal form”. The analysis of the dynamics of the normal forms 
reveals a qualitative picture of the flows of each bifurcation type. 
For the sake of clarity, Eq. (36) is written again as below: 

v c = aiv c -a>- l v c +f(v c , v c ) ^ 

v c = a^v c + oti v c +g(v c , v c ) 

in which with ~\v c is the complex conjugate of v c . 

As set of Eq. (37) are a pair of complex conjugate equations; one 
single equation needs to be studied eventually. A Flopf bifurcation 
is identified here and the result of these successive transforms is 
known [ 12 ]. 

v c = u A v c - ft)! v c + (C, v c - C 2 v c )(v 2 + v 2 c ) (38) 

with: 

16Ci = (f xxx "T fxyy + g X xy + Syyy) + l/ tWl \fxy (fxx + fyy) ~ gxyigxx 
"T Syy ) ~ fxxgxx fyygyy] 


16C 2 


(fxxy + fyyy ~ g xyy ~ g xxx ) + 1/^1 [2 (fxx +/xy) + 5(fyy 


+ 5g*x) + 5/ xx (fyy - g xy ) + 2 g 2 xy 

2gyy ~ fxy (gxx ^g>yy) 1 


fyygxy "f" ^gxxg 


yy 



In which f x (g x ) is the partial derivative to the first order of the func¬ 
tion /(g) with respect to the variable x and f xx (g xx ) is the partial 
derivative to the second order of a function / (g) with respect to 
the variable x. Using polar form ( v c = rcos 6 and v c = r sin 9, we final¬ 
ly obtain: 

r = air + Cir 3 

, (41) 

9 = coi r + C 2 r 


The amplitude of a limit cycle can be assed by the first of the 
previous equation. The non trivial solution establishes the exis¬ 
tence condition for the limit cycle as: ai/Ci < 0. The stability of 
each fixed point can be studied by the sign of the Jacobian J(r). 
From (41) J(r) = oc-i + 3C^ t therefore J(r fl ) = a 1 and /(r^) =-2ai. 
Thus, if ol\ is lower than zero and a greater than zero, the only sta¬ 
ble fixed point is r f i = 0 and if oq > 0 and a < 0 , the unique stable 
fixed point is and a limit cycle occurs. 


x p and x p [mm/s] 



x d basis for thermo-analysis - x p basis for thermo-analysis 

x d numerical results from dyn-analysis - x p numerical results from dyn-analysis 

x d semi-analytical results from dyn-analysis - x p semi-analytical from dyn-analysis 



Fig. 11. Piston and displacer movements. 
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Fig. 12. plots the iterative strategy used for the global analysis. 


Fig. 11 plots the results of the Eq. (27) using a direct numerical 
solution and the semi analytical approach previously described. As 
a comparison, the displacement of the piston and displacer from 
RE-1000 experimental data #1011 are plotted. 

Results for numerical and semi-analytical approaches match. 
The strokes are well determined (see Fig. lib) but the phase angle 
Uxdxp and the frequency show slight discrepancy. Comparisons de¬ 
tails will be given in next section. 

3. Model validation-RE-1000 comparison 

In order to validate the proposed modelling strategy, thor¬ 
oughly documented FPSE is needed. One of the available docu¬ 
mented FPSE is the RE-1000 studied by the NASA during the 80s 
[16]. The purpose of the tests was to collect data over a wide range 
of operating conditions. Changes in the operating condition were 
accomplished by varying the working fluid, the mean pressure le¬ 
vel, the heat source and sink temperatures. A dashpot load enabled 
to define the piston amplitude and displacer, pressure amplitudes 
and the engine frequency eventually. The description of one of the 
selected test configuration is given hereafter in Table 1. Some of 
the relevant variables (e.g. V sw , a, K) for these run are deduced 
from the raw data with Eqs. (3)-(5). 

3.1. Iterative evaluation strategy 

Thermodynamic and dynamic results are obtained using a two 
steps iterative strategy. Fig. 12 describes the general iterative 
scheme. The strategy used to couple the two models is first to ob¬ 
tain an evaluation of the chambers temperatures thanks to ther¬ 
modynamic analysis using kinematic parameters from RE-1000 


data. Second, these values are input parameters for dynamic anal¬ 
ysis for which the linear viscous damping coefficients d d , d p and 
dioad are used as calibration variables. They are modified to match 
the strokes of the piston and displacer as well as their phase angle 
in addition to the operating frequency. As defined by Eqs. (3)-(5), 
swept volume, phase angle and swept volume ratio can be evalu¬ 
ated eventually. 

After a first iteration for the thermodynamic plus dynamic anal¬ 
ysis, the evaluated kinematic results are parameters for a next iter¬ 
ation of thermodynamic analysis. The resulting temperatures T E 
and T c are compared with the previous values. A relative difference 
less than s T = 0.5% ensures the convergence of the thermodynamic 
analysis. This convergence criterion is denoted 7Y I+1 - 7V/ ^ s T in 
Fig. 13. 

The Reynolds values and the deduced pressure drop coefficients 
are strongly dependent on the kinematic variables as can be seen 
from (20). As a consequence, a second dynamic analysis is per¬ 
formed using the new thermodynamic results as parameters. 
New deduced operating values (V sw , k , a and co) are evaluated 
using Eqs. (3)-(5). These kinematic variables are defined as I(V. A 
convergence criterion is defined as: I<V i+ 1 - I< vi ^ s T . The maximal 
difference between successive results have to be less than s K = 2%. 

It is worthy of note that the pressure drops are not considered 
here as calibration parameters. Moreover, save for the very first 
iteration of all the comparison cases, the piston load is the single 
dynamic calibration parameter. 

3.2. RE-WOO run #1011 comparison 

As a first example, a detailed presentation is given for experi¬ 
mental run #1011. 
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Fig. 13. Clapeyron diagrams for compression and expansion spaces, (a) first iteration, and (b) second iteration. 


Table 2 

Results from thermodynamic analysis i = 0. 


Heat source T H (°C) 

600 

Cold source T c (°C) 

25 

Expansion temperature 

554.4 

Compression 

66.3 

T e (°C) 

(-0.1%) 

temperature T c (°C) 

(1.6%) 

Phase angle a (°) 

106.4 

Temperature ratio t (-) 

0.41 

(-1.7%) 

Swept volume V sw (cm 3 ) 

74.90 

Swept volume ratio k (-) 

1.022 

Operating frequency 

30.1 

Mean pressure p mea n (MPa) 

7.03 


/(Hz) 


Table 3 

Results from dynamic analysis from first (i = 0) thermodynamic analysis results 
parameters. 


Heat source T H (°C) 

600 

Cold source T c (°C) 

25 

Expansion temperature 

560.8 

Compression 

66.3 

Te (°C) 


temperature T c (°C) 


Phase angle a (°) 

108.11 

(-1.6%) 

Temperature ratio t (-) 

0.41 

Swept volume 

74.95 

Swept volume ratio k (-) 

1.06 

Vsw (cm 3 ) 

(0.06%) 


(-4.4%) 

Operating frequency 

27.8 (7.7%) 

Mean pressure p mecin 

7.03 

/(Hz) 


(MPa) 



• First iteration: Thermodynamic analysis 

Input data and results from thermodynamic analysis first eval¬ 
uation are given in Table 2. Bold figures represent the estimated 
values from the model whereas input data are given in gray back¬ 
ground for information. For each calculated value, relative devia¬ 
tion with the corresponding experimental reference value is 
given in brackets. Temperature evaluations are in good agreement 
with experimental ones. 

• First iteration: Dynamic analysis 

Input data and results are given in Table 3. Bold figures repre¬ 
sent the estimated values from the model whereas input data from 
the thermodynamic results are given for information. 

The output power is evaluated through the linear damping of 
the piston: P outdyn = 1166 W. Compare to the RE-1000 output 
power Pout = 955 W the model gives a first quite good estimation 
(22% discrepancy). 

Fig. 13a) plots Clapeyron diagrams for the expansion and the 
compression chamber. The results from the dynamic and the ther¬ 


Table 4 

Results from thermodynamic analysis (i = 1). 


Heat source T H (°C) 

600 

Cold source T c (°C) 

25 

Expansion 

552.9 

Compression 

65.5 

temperature T E (° C) 

(0.18%) 

temperature T c (°C) 

(0.26%) 

Phase angle a (°) 

108.11 

Temperature ratio t (-) 

0.41 




(0.08%) 

Swept volume V sw (cm 3 ) 

74.09 

Swept volume ratio k (-) 

1.06 

Operating frequency 

27.8 

Mean pressure p mean 

7.03 

/(Hz) 


(MPa) 



Table 5 

Results from dynamic analysis (/=!). 


Heat source T H (°C) 

600 

Cold source T c (°C) 

25 

Expansion 

552.9 

Compression 

65.5 

temperature T E (°C) 
Phase angle a (°) 

110.3 

temperature T c (°C) 
Temperature ratio t (-) 

0.41 

Swept volume V sw 

(-2.0%) 

72.60 

Swept volume ratio 

1.08 

(cm 3 ) 

(1.7%) 

K (-) 

(-1.8%) 

Operating frequency 

27.8 

Mean pressure p mean 

7.03 

/(Hz) 

(-0.3%) 

(MPa) 



modynamic first analysis are presented and appear to be very close 
to each other. 

Though the convergence criteria is satisfied, new values for 
Reynolds number have to be evaluated again in order to take into 
account drift related to swept volume ratio and operating fre¬ 
quency. Modified pressure drops values have strong influence on 
the nonlinear characteristics of the system. As a consequence the 
dissipative load variable will be adapted and so do the output 
power. 

• Second iteration: Thermodynamic analysis 

For the second iteration, kinematic results from the dynamic 
analysis are set as input parameters. Bold figures of Table 4 rep¬ 
resent the estimated values from the model whereas input data 
from the previous iteration are given for information. For each 
calculated value, the drift with the previous evaluation is also 
given. 

Temperature drifts are small which ensures convergence 

(e T <0.5%). 
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Table 6 

Final comparison for run #1011. 


Expansion temperature T E (°C) 

553.2 (-0.27%) 

Compression temperature T c (°C) 

65.0 (1.17%) 

Piston stroke (cm) 

2.72 (-2.6%) 

Displacer stroke (cm) 

2.41 (-2%) 

Displacer - piston phase angle (°) 

62.2 (8.2%) 



Phase angle a (°) 

110.3 (3.64%) 



Swept volume V sw (cm 3 ) 

72.60 (-3.08%) 

Swept volume ratio k (-) 

1.08 (6.16%) 

Operating frequency/(Hz) 

27.9 (-7.47%) 

Output power P out (kW) 

1.1 (15.60%) 

Overall efficiency rj tm (%) 

29.2 (23.3%) 

Input thermal power * (kW) 

3777 (-6.46%) 



Output thermal power (kW) 

2107 (-30.5%) 


* For the experimental results, the electrical power input is used as an evaluation of the power input. 


Table 7 

Results from different models of RE-1000. 


Parameter 

Urieli and 
Berchowitz. 

Walker 

and 

Senft 

Rogdakis 

case 

study 1 

Rogdakis 

case 

study 2 

Proposed 

global 

analysis 

Operating frequency 

10.7 

9.7 

0 

0.3 

-7.5 

/(%) 

Displacer - piston 

36.2 

29.6 

28.9 

29.4 

8.2 

phase angle (%) 

Amplitude ratio (%) 

-41.5 

-40.6 

-17.9 

-31.1 

-0.9 


• Second iteration: Dynamic analysis 

From the results of from the dynamic analysis given in Table 5, 
the output power evaluated through the linear damping of the pis¬ 
ton: P out = 1104 W. Compare to the RE-1000 output power P ou texp = 
955 W the model gives a good estimation (15.6% discrepancy). 

It is worthy of note that the dissipated power at the displacer is 
predicted as about 114 W which is about 10% of the piston power. 
As discrepancies for kinematic variables show less than 2% differ¬ 
ence with the first iteration the coupling of the approaches appears 
to be a stable process. After the second iteration, Clapeyron dia¬ 
grams for the expansion and compression chambers show a good 
agreement in Fig. 13b). 

• Final comparison with run #1011 data 

Table 6 gives the final results after the two steps procedure for 
the global analysis. 

The proposed global analysis allows an estimation of the kine¬ 
matic results close to 10%. Compare to linear dynamic analysis, 
the global approach gives better results and is able to evaluate 
thermodynamic as well as dynamic variables. 

Table 7 gives comparisons of dynamical models from [6] with 
experimental results. For the chosen kinematic results, the 


Pout [kW] 



Table 8 

Parameter values of run #1360. 


Heat source T H (°C) 

600 

Cold source T c (°C) 

25 

Phase angle a (°) 

111.057 

Swept volume ratio k (-) 

1.605 

Swept volume V sw (cm 3 ) 

83.7476 

Working fluid 

Helium 

Mean pressure p mea n (MPa) 

7.038 

Overall efficiency r/ tm (%) 

21.4 

Operating frequency/(Hz) 

31.1 

Output power (W) 

987 

Wire mesh regenerator 1 


regenerator length (mm) 

64.46 

Wire diameter d w (pm) 

88.9 

Porosity of matrix % (-) 

0.759 

Regenerator hydraulic 

0.07 

Material conductivity 

16.3 

radius r hR (mm) 


(Witt 1 K _1 ) 


Displacer 2 


Rod diameter (mm) 

18.085 

Displacer weight (g) 

381 

Gas spring volume (cm 3 ) 

18.8 


proposed model appears to be closer than any of the approaches. 
Because the evaluated temperatures are also close to the experi¬ 
mental ones (see Table 6), performances of the engine can be well 
predicted eventually. 

3.3. Influence of piston stroke amplitude for displacer 1 

Experimental procedure used in [ 16] consists in defining the pis¬ 
ton stroke. From run #1006 to #1012, piston stroke varies from 2 to 
3 cm with 0.2 cm increment by load adjustment. By increasing the 
piston stroke, the displacer stroke as well as the phase angle are dif¬ 
ferent for each test. Finally, swept volumes increase with constant 
operating frequencies. The global model is used to determine the 
corresponding responses. Fig. 14 plots the output power and overall 
efficiency with respect to the swept volumes. The evolutions of the 
performances are well represented with a 14.1 % mean departure for 
powers, 21.8% for efficiencies and 3.2% for swept volumes. Standard 
deviations are 2.4, 2.4 and 1.1 respectively. 

3.4. Influence of piston stroke amplitude for displacer 2 

In addition to the evaluation of the effect of the piston stroke 
amplitude for one configuration (runs #1006 to #1012) another 


n [%] 



55 60 65 70 75 80 

(b) 


Fig. 14. (a) output power, and (b) overall efficiency as a function of the swept volume for displacer 1. 
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Pout [kW] 



(a) 



Fig. 15. (a) output power, and (b) overall efficiency as a function of the swept volume for displacer 2. 


configuration have been chosen (run #1353 to #1359) for model 
validation. For these tests, the displacer mass is reduced to 381 g. 
Table 8 gives some of the parameters. It is worthy of note that 
the length of the lighter displacer 2 is less than displacer 1. By 
doing this, dead volumes of the new engine configuration are more 
important. 

The evolutions of the performances are well represented as can 
be seen in Fig. 15. The mean departure for powers is 27.9%, 36.6% 
for efficiencies and 0.8% for swept volumes. Standard deviations 
are 7.3, 2.5 and 3.6 respectively. 

4. Conclusion 

Semi analytical thermodynamic and dynamic models of FPSE 
have been elaborated. The models integrate the regenerator effec¬ 
tiveness, heat exchangers performances as well as conduction 
losses. Besides a suitable analysis of the pressure drop losses is 
conducted since the steady state characteristics strongly depends 
on these effects. Both thermodynamic and dynamic evaluation 
can be performed in the developed global analysis strategy. The 
available RE-1000 experimental data have been compared to the 
model results which show good agreements. Because the geomet¬ 
rical, thermal, fluid and dynamic properties are direct parameters 
of the global approach, it can be used as a preliminary design tool 
for FPSE or guidelines for optimization of existing Stirling engines. 
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